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As atoms formed for the first time during primordial recombination, they emitted bound-bound 
and free-bound radiation leading to spectral distortions to the cosmic microwave background. These 
distortions might become observable in the future with high-sensitivity spectrometers, and pro- 
vide a new window into physical conditions in the early universe. The standard multilevel atom 
method habitually used to compute the recombination spectrum is computationally expensive, im- 
peding a detailed quantitative exploration of the information contained in spectral distortions thus 
far. In this work it is shown that the emissivity in optically thin allowed transitions can be fac- 
tored into a computationally expensive but cosmology-independent part and a computationally 
cheap, cosmology-dependent part. The slow part of the computation consists in pre-computing 
temperature-dependent effective "conductances", linearly relating line or continuum intensity to 
departures from Saha equilibrium of the lowest-order excited states (2s and 2p), that can be seen as 
"voltages" . The computation of these departures from equilibrium as a function of redshift is itself 
very fast, thanks to the effective multilevel atom method introduced in an earlier work. With this 
factorization, the recurring cost of a single computation of the recombination spectrum is only a 
fraction of a second on a standard laptop, more than four orders of magnitude shorter than standard 
computations. The spectrum from helium recombination can be efficiently computed in an identical 
way, and a fast code computing the full primordial recombination spectrum with this method will 
be made publicly available soon. 



I. INTRODUCTION 

A significant part of our knowledge about the universe 
at early times and on large distance scales is derived 
from the observation and analysis of its spatial inhomo- 
geneities. In particular, observations of the temperature 
and polarization anisotropies of the cosmic microwave 
background (CMB) have allowed cosmologists to deter- 
mine the geometry, contents, and initial conditions of the 
universe to an exquisite level of precision (see for example 
Ref. pQ). 

Additional, and perhaps complementary information, 
is hidden in the tiny but unavoidable spectral distortions 
to the nearly perfectly thermal radiation background 
[21 [3]. On the one hand, broad spectral distortions can be 
generated by energy injection in the early universe, tak- 
ing the form of chemical potential (fj,— type) or Compton 
(y~ type) distortions (and more generally continuously 
interpolating between these two analytic cases [4]), de- 
pending on the redshift of energy injection. Physical 
mechanisms causing such energy injections include the 
dissipation of small-scale acoustic waves, which occurs in 
the standard cosmological picture (see e.g. [5]), and pos- 
sibly the decay or annihilation of dark matter into stan- 
dard model particles that then deposit their energy into 
the primeval plasma [6]. A vast literature exists on the 
subject, and we refer the interested reader to the recent 
works [H [3 |H] and references therein. 

On the other hand, the primordial recombinations of 
helium and hydrogen lead to a few distortion photons 
per atom, in the form of free-bound and bound-bound 



photons. The seminal papers on primordial recombina- 
tion by Peebles and Zeldovich et al. [9l [10] already eval- 
uated the distortion due to Ly-a transitions and 2s — Is 
two-photon decays. Even though this represents a large 
distortion to the Wien tail of the CMB blackbody spec- 
trum, it lies many orders of magnitude below the cosmic 
infrared background (CIB, [H]), which renders its detec- 
tion very challenging, if not hopeless. 

In addition, exactly one free-bound photon per hydro- 
gen atom and two per helium atom were emitted, as well 
as a few bound-bound photons per atom from transi- 
tions between highly excited state^ [12]. It was first 
pointed out by Dubrovich [T3] that these primordial re- 
combination lines could be observed today as broadened 
features in the centimeter to decimeter wavelength range. 
A few authors have since then tackled the computation 
of the recombination spectrum, with various degrees of 
approximation or numerical convergence |14H17j . and it 
is only quite recently that highly-accurate computations 
were carried out by Rubino-Martm, Chluba and Sunyaev 



Thus far the study of the primordial recombination 
spectrum has remained the niche domain of a few afi- 
cionados. First and foremost, the observational prospects 
may seem meager, as the predicted signal in the cm- 
wavelengths lies a billion times below the undistorted 
CMB spectrum, far below the current best upper bounds 
on spectral distortions from FIRAS [21] , In addition, 
galactic and extragalactic foregrounds would need to be 
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1 Throughout this paper we refer to bound-bound transitions be- 
tween "highly excited" states as those for which the lower state 
itself is excited, n' — > n > 2. 
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understood and subtracted very precisely. Finally, the 
machinery required to compute a recombination spec- 
trum with the standard multilevel method is, although 
conceptually not very difficult, somewhat cumbersome to 
implement and too computationally expensive for a sys- 
tematic analysis of its information content. 

The reward in finding such a needle-in-a-haystack sig- 
nal is, however, potentially significant. Ref. [22] showed 
that the recombination spectrum is a sensitive thermome- 
ter and baryometer. It could also provide a clean mea- 
surement of the primordial helium abundance, before the 
formation of the first stars [3J [20] . Finally, pre-existing 
spectral distortions could lead to a significant increase 
of the recombination radiation even if the initial distor- 
tions are small in absolute value [T5J [23J . The recombina- 
tion spectrum could therefore be a probe of non-standard 
physics such as dark-matter annihilations [24] . 

Let us point out, in addition, that technological ad- 
vances should make it possible to reach sensitivities three 
orders of magnitude below that of FIRAS, corresponding 
to distortions at the level of ~ 1CP 8 (see the proposal for 
the PIXIE instrument [2"5]). Spectral distortions from re- 
combination are only one order of magnitude weaker than 
this sensitivity (and even get to the 10 -8 level around 10 
GHz [12]), and it is not unlikely that they will be within 
reach of the next generation instruments. The feasibility 
of foreground subtraction at the required level has yet 
to be demonstrated, but one may hope that the spatial 
isotropy of the signal and its very specific spectral fea- 
tures should allow us to disentangle it from foreground 
emission. 

Before any of the truly challenging issues of instrumen- 
tal sensitivity and foreground subtraction are addressed, 
it seems that the first task is to undertake a detailed 
quantitative study of the information content of the re- 
combination spectrum. In order to do so, a fast and 
accurate computational method is required, so that the 
large space of cosmological parameters can be efficiently 
explored^ Introducing such a fast method is the purpose 
of the present work. 

The main task in obtaining the spectrum resides in 
computing the populations of the excited states, or, more 
precisely, their small departures from equilibrium with 
one another, since the line and continuum emissions are 
proportional to the latter. High-precision spectra require 
accounting for excited states up to principal quantum 
number n max of a few hundred, resolving the angular mo- 
mentum substates. With the standard multilevel atom 
method, one has to invert a large "° ax x matrix 
at each timestep in order to compute the populations of 
the excited states (although this matrix is sparse due to 
selection rules, so only C(n^ ax ) elements are nonzero). 



2 It was brought to my attention by J. Chluba that Fendt I2GI con- 
ducted a preliminary study of cosmological parameter estimation 
from spectral distortions, using the fast interpolation algorithm 
Pico [27]. 



Currently the fastest code using this method takes about 
one hour for n max = 100 and one full day for n max = 250 
on a standard laptop, with the computation time scaling 
as tocn&L M- 

The method that we introduce here allows us to fac- 
torize the problem into a computationally expensive part 
(for which large linear systems need to be solved) that 
is cosmology-independent and can be pre-computed once 
and for all, and a computationally cheap part that does 
depend on cosmology. This method builds on and ex- 
tends the effective multilevel atom (EMLA) method in- 
troduced in an earlier work [2H] (hereafter AH10; see also 
Refs. 30 Siy, it is, however, not a trivial extension, 
since the EMLA method is designed for computing the 
free electron fraction x e (z) and essentially collapses all 
the transitions between excited states into effective tran- 
sition rates into and out of 2s and 2p. In this process, 
information not directly necessary to the evolution of the 
free electron fraction is lost, whereas our present goal is 
to go beyond x e (z) and compute the full recombination 
spectrum. 

We shall lay down our method in detail in the remain- 
der of this paper, but the main idea can be intuitively 
understood if one pictures the system of radiatively con- 
nected excited levels as a circuit, where the currents are 
the line intensities, the voltages are the departures of the 
excited states populations from Saha equilibrium, and 
the conductances are the transition rates (this insight- 
ful analogy is due to Chris Hirata [32]). The linearity of 
Kirchhoff 's laws (the steady-state rate equations for the 
excited state) ensures that the "current" in any transi- 
tion is proportional to the outer "voltages" , i.e. the de- 
partures from Saha equilibrium of the 2s and 2p states. 
The proportionality coefficients, which we shall call effec- 
tive conductances, moreover only depend on the temper- 
ature of the ambient blackbody radiation that is nearly 
undistorted at the relevant frequencies (this can in princi- 
ple be generalized to include simple parameterizations of 
the ambient spectrum, as well as collisional transitions). 
Once the effective conductances are pre-computed as a 
function of transition energy and temperature, the re- 
combination spectrum can be computed very efficiently 
for any cosmology, by using the EMLA method to eval- 
uate the recombination history and the outer "voltages" 
as a function of redshift. This method is very similar in 
spirit to the widely used line-of-sight integration method 
for CMB anisotropics [33], which factorizes the compu- 
tation of the CMB power-spectrum into a geometric, 
cosmology-independent part and a cosmology-dependent 
but multipole-independent source term. We illustrate our 
method graphically in Fig. [I] 

This paper is organized as follows. In Sec. [Hj we write 
down the general equations that need to be solved for the 
computation of the recombination spectrum. The effec- 
tive conductance method is described in Sec. Mil We dis- 
cuss our numerical implementation and results in Sec. |IV| 
and give our conclusions and future research directions 
in Sec. El 
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FIG. 1. Schematic representation of the main idea of this 
paper, highlighting the circuit analogy of Hirata [32]. All 
transitions between excited states are mediated by blackbody 
photons and their rates (and corresponding resistances R(T r ) 
shown as solid resistor symbols) only depend on the blackbody 
temperature. The linearity of the steady-state rate equations 
(Kirchhoff 's laws in the circuit analogy) ensures that the cur- 
rent I n ' n in any transition between excited states is linear in 
the departures from Saha equilibrium Ax2s and Ax2 P , which 
play the role of externally imposed voltages. The proportion- 
ality coefficient is a function of temperature only. The values 
of Ax2s and Ax2 P are determined from the balance between 
the net downward current from the highly excited states (ob- 
tained through the effective rates defined in AH10) and the 
non-thermal transition rates to the ground state (represented 
by dashed resistor symbols). We have not represented all al- 
lowed transitions in order to keep the graph readable. 



II. GENERAL EQUATIONS 



A. Notation 



We denote by «h the total number density of hydro- 
gen in all its forms (ionized and neutral), x e the ratio of 
the free-electron abundance to the total hydrogen abun- 
dance, x p the fraction of ionized hydrogen, and x n i (or 
in some cases X n {) the fractional abundance of neutral 
hydrogen in the excited state of principal quantum num- 
ber n and angular momentum number I. The matter 
and radiation temperatures are denoted by T m and T x , 
respectively. We denote emissivities by j v (with units 
of energy per unit time per frequency interval per unit 
volume per unit solid angle) and specific intensities by I v 
(with units of energy per unit time per frequency interval 
per unit area per unit solid angle). All our derivations 
are for hydrogen atoms but the generalization to helium 
is straightforward. We only consider radiative transitions 
here and neglect the effect of collisions. 



B. Bound-bound emission from transitions 
between excited states 

The emissivity due to bound-bound transitions be- 
tween excited states is given by 

]hh{v) = "H^j— X 

^ ^ [Xn'l'Rn'V^Hd - XnlRnl^n'l'] S{v - *Vn)i(l) 
2<n<n' 1,1' 

where v n i n is the frequency of the n'V — > nl transitions 
and Rni-yn'V represents the radiative transition rate from 
nl to n'V . The recombination process adds at most a few 
photons per atom, and the transitions between excited 
states are mostly below the peak of the blackbody spec- 
trum (except for Balmer transitions, but their energy is 
only a few times above the blackbody peak, where there 
is still a very large number of thermal photons per hydro- 
gen atom). As a consequence, the radiation field medi- 
ating the transitions is, to ~ 10~ 8 accuracy, a blackbody 
at temperature T T . One can check the validity of this as- 
sumption a posteriori once the distortions are computed 
(see for example Fig. 2 of Ref. [H]). Note, however, that 
small y- distort ions to the blackbody spectrum (at the 
level of y ~ 10~ 6 ) may significantly enhance the hydro- 
gen and helium line emission [TSJ [23] ■ We defer the study 
of the effect of such pre-existing distortions on the recom- 
bination spectrum to future work, and shall here assume 
that transitions between excited states are mediated by 
thermal photons only. 

With this assumption, the radiative transition rates 
satisfy the detailed balance relations, 

Rn'l'^nl = — —Rnl-yn'l', (2) 
Qn'l' 

where we have defined 

g„ J = (2l + l)e- B -/ fc3i , (3) 

where E n = —13.6 eV/n 2 is the (negative) energy of the 
bound states with principal quantum number n. 

We see that if the excited states were in Boltzmann 
equilibrium, so that x n >i>/x n i = q n 'i'/<lnli the net emis- 
sivity would vanish. The emissivity therefore scales lin- 
early with the small departures from equilibrium. To 
make this apparent, let us define the departures from 
Saha equilibrium at temperature T r with the free elec- 
tron and protons, 

Ax n i = x n i - —n H x e x p , (4) 

Qe 

where we have defined 

/27rm e fcT r \ 3/2 
q e = [—f^- ) , (5) 

where rn e is the reduced mass of the electron-proton sys- 
tem. We can now rewrite the bound-bound emissivity in 
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the following form: 

Jbb(^) = n n ~. x 
An 



E E 



q n i 
q n 'l> 



Ax r 



•v 



Ax, 



2<n<n' 1,1' 

where the Saha-equilibrium pieces have cancelled out. 
C. Free-bound emission 



free-bound emissivity would vanish if the excited states 
were in Saha equilibrium with the ionized plasma and if 
the matter and radiation temperature were identical. We 
can rewrite the free-bound emissivity in terms of small 
departures from equilibrium as follows: 



jfbO) = 

EE 

n>2 Kn 



hv 

qm , -.AT 

— n,Yix e x p -i n {v) — 

Qe 1 



Ax n l 



dv 



(10) 



The emissivity due to free-bound transitions to excited 
states is given by 



hv 



Jfb(") = nH^X)X! 



n>2 Kn 



da 



riY{X e x 



ill 



p dv 



dp n l 

dp 



> (7) 



where da n i/dv is the differential recombination coeffi- 
cient per frequency interval of the emitted photon and 
dfi n i/dv is the differential photoionization rate per fre- 
quency interval of the ionizing photon. Recombinations 
of the thermal electrons and protons to the excited states 
are mediated by blackbody photons; as a consequence, 
da n i/dv and dj3 n i/dv are related through the relation 
(see for example Eq. (2) of AH10): 



da 



Qnl d(3 n i 
Q P dv 



3/2 



exp 



{hv + E n 



1 



1 



kT-n 



(8) 



We have purposefully made the ratio and difference of 
the matter and radiation temperatures appear in this 
expression. At high redshifts when the recombination 
spectrum is emitted, the matter temperature is locked 
to the radiation temperature by Compton heating (see 
for example Ref. [34]). Computing the coupled evolu- 
tion of the ionization history and matter temperature 
with HyRec^] [35], we find that the fractional differ- 
ence between the two temperatures is less than 10~ 5 for 
z > 1100 and less than 10~ 3 for z > 800. We can there- 
fore Taylor-expand the above expression in the small pa- 
rameter AT/T = (1 - T m /T T ) (note that with this con- 
vention AT > for T m < T r ) and obtain 



da nt q„i d(3 n i 



dv 



q e dv 

Qnl d(3 n l 



1 



h{v 



kT T 



AT 



Qe 



dv 



AT 



(9) 



where v cn = —E n /h is the photoionization threshold 
from the n-th shell and the second expression defines the 
dimensionless parameter j n (v). Again, we see that the 
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D. Radiative transfer equation 

Strictly speaking, the computation of the specific in- 
tensity requires the knowledge of not only the emissivity 
but also the absorption coefficient [36]. However, the 
Sobolev optical depths of most bound-bound transitions 
are much lower than unity |371 I38j . perhaps with the 
exception of the very high-rt transitions (n > 300), as 
can be seen from extrapolating Fig. 11 of Ref. [19]. It 
is possible that a proper accounting of the nonzero opti- 
cal depth in very high-n transitions could lead to small 
modifications in the low-frequency part of the spectrum 
(y < 0.1 GHz); however, in that frequency range other 
effects that we are neglecting significantly affect the spec- 
trum too, such as free-free absorption [19] and collisional 
transitions [21] . We shall therefore assume the optically- 
thin limit for all bound-bound transitions between ex- 
cited states, as well as free-bound transitions. 

The radiative transfer equation in the optically thin 
regime in an expanding homogeneous universe takes the 
simple form 



d 

dt 



Iu 



d_ 
di 



Hv 



d_ fl, 
dv 



3v 



(11) 



In between resonances (where j v = 0), the quantity I^/v 3 
is conserved along a photon trajectory, so that 

I v {z)=(^f (I 2 ) 



where 



l + z' = -(1 + z). 



(13) 



In the vicinity of a resonance line, j u — JqS(v — v ), the 
solution to the radiative transfer equation is 

cJ 



12 



Hv -' 



(14) 



where Ij and I~ are the specific intensities at the blue 
and red sides of the resonant line, respectively. 

The general solution of the radiative transfer equation 
can be written in the following integral form 



I v (z) = c I dv' 

J v 



l + z 
l + z' 
~ J dv' Jv , 



Hv 1 



Hv' njj 



(15) 
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where the redshift z' and frequency v' are related through 
Eq. (13) and we used the fact that njj(z) oc (1 + z) 3 . 



III. THE EFFECTIVE CONDUCTANCE 
METHOD 

A. Populations of the excited states 

In order to obtain the bound-bound and free-bound 
cmissivities, we see that we need to evaluate the popula- 
tions of the excited states, or, more precisely, their small 
departures from Saha equilibrium with the plasma. 

The populations of excited states can be obtained to 
very high accuracy in the steady-state approximation, be- 
cause of the large ratio of internal transition rates to the 
overall recombination rate. This assumption was checked 
explicitly in Ref. [25] and found to be extremely accurate, 
for the computation of both the recombination history 
and the recombination spectrum. 

Following AH10, we separate the excited states in "in- 
terior" states, only connected radiatively to other excited 
states and the continuum, and "interface" states, essen- 
tially 2s and 2p (and potentially any additional "weak 
interface" states, such as 3s, 3p, 3d...), which are radia- 
tively connected to the ground state. We denote the pop- 
ulations of the former by a capital Xk, where K stands 
for both quantum numbers of the state, and those of the 
latter by Xi, where i = 2s, 2p (and 3s, 3p, 3d... if needed). 

The steady-state rate equation for the population of 
the interior state K is 

M X K = nftX e X p aK + ^ XiRi^K 

i 

+ x lRl^k - X K T Kl (16) 

where (T m , T r ) is the total recombination coefficient to 
the state K (accounting for stimulated recombinations) 
and 



K 



(17) 



is the total rate of transitions out of the state K, where 
is the total photoionization rate from K. Note that 
here again we have assumed that the Sobolev escape 
probability is unity in all transitions (or that the Sobolev 
optical depth is zero). If this were not the case the net 
transition rates would depend non-linearly on the state 
populations, which would significantly complicate mat- 
ters. 

Provided the transitions between excited states are me- 
diated by blackbody photons, the transition rates satisfy 
the detailed balance relation Eq. The recombination 
coefficients and photoionization rates are related through 



<le 



(18) 



We can now rewrite the system in terms of the small 
departures from Saha equilibrium with the continuum, 
and to linear order in AT/T: 







-ri}ix e XpAT 



da 



K 



T m =T r 



AX L R L ^ K - AX K T 



K ■ 



(19) 



In the standard multilevel atom method, the large sys- 
tem ( 19 ) is solved at every timestep, which makes the 



computation very slow. 

Following AH10, we define the matrix M of elements 



M KL = T K 8 KL - (1 - 5 KL )R K ^ L . 



(20) 



We also define the vector AX of elements AXk and the 
source vector AS of elements 



ASk = —nnx e XpAT 



da K 



T m =T r 



£ Ax % R t 



(21) 



The system ( 19 1 can be rewritten in compact matrix form 
as 



M T (AX) = AS, 



(22) 



where M T is the transpose of M. We showed in AH10 
that the matrix M is nonsingular, and this system has 
the formal solution 

AX = (M T )- 1 (AS') = (M- 1 ) T (AS'), (23) 

i.e., explicitly, 

AX K = ^(M-^lkASl 

L 

= -Ti ii x e x p ATy^(M~ 1 ) LK '^- 

dT m T m =T r 

+ ^Ax l ^(M- 1 ) LK iW. (24) 

t L 

We showed in AH10 that detailed balance relations be- 
tween radiative transition rates ensure that 



{M~ l ) LK = —(M~ 1 )kl- 



(25) 



Using the detailed balance relation for R^l, we may 
rewrite the last term of Eq. (124]) as 



Qi r 



- IK p . 
= r\ 



K i 



(26) 



where P\ is the probability that an excited atom initially 
in the interior state K eventually reaches the interface 
state i (after possibly many transitions in the interior), 
the complementary events being to eventually reach one 
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of the other interface states or to be photoionized. The 
probabilities P % K were defined in AH10, where they were 
an intermediate step to compute the effective recombina- 
tion coefficients to and effective transition rates between 
the interface states. In the last line of Eq. (26 1, we have 



used the formal solution of the defining equation for the 
probabilities PL. 

Let us now deal with the first term of Eq. ( 24 1 . We 
define the dimensionless coefficients 



Ik 



d log a K 



d\ogT u 



(27) 



Using again detailed balance relations, we obtain 



da i 



}LK 



T m =T r 



L 



T r q e 



1 fen 



P, 



K - 



1 )klPl1l 



(28) 



where the last equality defines the dimensionless coeffi- 
cient Pjj-. If all the coefficients 7^ were equal to unity, 
then we would have P|- = P^, the probability that an 
excited atom initially in the interior state K eventually 
gets photoionized before reaching an interface state. In 
general, however, 7l 7^ I (but is in general positive and 
of order unity), so the numbers Pjj do not have a clear 
physical significance but are numerically of the same or- 
der as the Pjr. 

We therefore end up with the following compact ex- 
pression for the departures from Saha equilibrium: 



where we have defined the coefficients 
q n i 



1,1' 



pe pe 

r n'V r nl 



Pi 



(32) 



Qi, n {T T ) [Kl ~ Pi>t>] Rn^n'V- (33) 

LI' qi 



31 ): we have chosen this con- 



Note the minus sign in Eq. 
vention because the excited states are in general under- 
populated with respect to Saha equilibrium, and the co- 
efficients Qn'n defined in Eq. (33) are positive (in general 
the probability of reaching interface states decreases as 
n increases). 

If one thinks of the departures from Saha equilibrium 
Axi as voltages and of the net decay rate in the n' —> n 
transitions as a current, then the coefficients Q n , n can be 
thought of as effective conductances linearly relating the 
two. A similar analogy can be made for the first term: 
there, the voltage is the fractional temperature difference 
AT/T and the conductance would be n\\x f ,XpQ^ l i n . Note 
that the Q e have units of cm 3 s _1 whereas the Q % have 
units of s _1 . 

Similarly, the net free-bound decay rate per unit fre- 
quency can be rewritten as 



EE 

n>2 Kn 



q n i , ,AT 

— nnx e x p ^ n {v)— Ax n i 

qe 



dv 



dgjj, AT 



dv 



AY 1 K 6e AT 

AX K = — P K nnx e x p —- 
q P T 



E 



P l K Ax t 



(29) where we have defined the differential effective conduc- 
tances 



Looking more closely at equation ( 29 1 , we see that the 



departure from Saha equilibrium of any excited state is 
proportional to the difference between matter and radia- 
tion temperature (so long as this difference is small) times 
nnx e x p , and to the departures from Saha equilibrium of 
the small set of interface states. The proportionality co- 
efficients are functions of radiation temperature only. If 
we define Pj = 5ij and Pf = for interface states we 
can write a general equation valid for any excited state 
(including the interface states) in the form 



dv 



dG 



EE? -*.(") 



n>2 Kn 



n>2 Kn 



qni n , dj3 n l 



dv 



df3, 



dv 



(35) 
(36) 



Here again we have defined the effective conductance 
dQl h /dv with a positive sign, such that the current is pos- 
itive when the interface states are underpopulated with 
respect to Saha equilibrium. 



a _ 1 nl w AT 

<-± x nl — r nl n\\X e Xp 

q& i 



\ - q n i 



q, 



PlrAxi 



(30) 



C. Bound-bound and free-bound emissivities 



B. Effective conductances 

We can now rewrite the net decay rate in the n' — > n 
transitions (with n < n'), per hydrogen atom, in the form 



E 



q n i 
q n 'v 



Ax n , v - AXnl 



AT 



Rnl—>n't 



= Gl^nnXeXp— - Gl'n Ax h (31) 



The expressions for the emissivities can now be rewrit- 
ten in terms of the effective conductances: 

. . , hv 



E 



AT 



g e n , nnil X e X p — - Y Gn'n A: > 



jfb(^) = «H 



47T 



dv 



AT 



n. ...... V-^ Ax , 



S(v - v n > n )(37) 
(.38) 



dv 
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We can rewrite the total emissivity formally as 

• / -> ■ f \ hv 
3v = Jbb(^) + 3fb\y) = n H-r- x 

47T 

dg e at ^dg* 

(v,T T ) n a x e x p — ) — — (i/,T r )Aa 



rZ;/ 



where 



e,i 
fb 



n<n' 



dv 



,(39) 



(40) 



Equation ( 39 1 , along with the definitions of effective con- 
ductances given above, constitutes the fundamental re- 
sult of this paper. What it means is that the emissiv- 
ity can be factored into a cosmology- independent part, 
embodied in the effective conductances, that is compu- 
tationally expensive but can be pretabulated as a func- 
tion of temperature, and a simple cosmology-dependent 
part, entering through the departures from Saha equilib- 
rium of 2s and 2p and temperature differences, which are 
straightforward to compute for each particular cosmology 
with the effective multilevel atom method developed in 
AH10. 

Before proceeding further, let us point out that even 
though we have Taylor-expanded our expressions in the 
small difference between matter and radiation temper- 
ature, this is not required. One can very easily obtain 
more general expressions for arbitrary AT/T, in which 
the coefficient dg e /dv would become a function of radia- 
tion and matter temperatures. We leave it as an exercice 
for the interested reader to derive the exact expressions. 



IV. IMPLEMENTATION AND RESULTS 
A. Computation of effective conductances 

The computation of the effective conductances requires 
solving large (formally infinite) linear systems. One must 
impose some truncation criterion to make the system fi- 
nite and tractable. Following previous works, we simply 
ignore all excited levels above some cutoff value of the 
principal quantum numbei|^]n max . 

We tabulated the effective conductances on a grid of 
temperatures for several values of n max ranging from 
60 to 500. Matrix elements for bound-bound transi- 
tions were computed as in Ref. 39 and those for free- 
bound transitions as in Ref. Hi 



Equations (261 and 



(281 for the probabilities P^-(T r ) and the dimension- 



of the problem, as they require solving large (of order 
n max/2 x n max/2) matrix equations. Due to selection 
rules for radiative transitions, these large matrices are 
very sparse, and only have of order n^ ax nonvanishing 
elements. We can therefore use a sparse matrix tech- 
nique identical to the one introduced in Ref. [38 . 

We show some of the computed effective conductances 
in Fig. [2j for a temperature T T = 3800 K (corresponding 
roughly to the peak of the emission) . 



B. Practical simplification 



Equation ( 39 1 is a rather remarkable result in its raw 



form, but its implementation without further simplifica- 
tion could be somewhat cumbersome: if considering ex- 
cited states up to principal quantum number n max , one 
would in principle need to tabulate of order n^ ax /2 func- 
tions g n 'n(T r ) as a function of temperature. For n max of 
a few hundred, required for a fully converged spectrum 
in the GHz region, one would need to store tens of thou- 
sands of coefficients on a fine grid of temperature values, 
interpolate them at each redshift, and compute the re- 
combination spectrum on a grid fine enough to resolve 
all the resonances. In addition, whereas hydrogen and 
singly ionized helium benefit from the accidental energy 
degeneracy between angular momentum substates, this is 
not the case for neutral helium, which as a consequence 
has a much larger set of lines. 

In order to save a significant amount of memory with 
negligible cost in accuracy (as we shall demonstrate in 
the next section), we group resonances into bins b of finite 
width Az/f,, so we use 



dg e >' 

dv 



where v\, are the bin centers, and 



n<n' 



dg 



fb 



dv 



(41) 



(42) 



less numbers P^(T T ) represent the time consuming part 



where lj(i') is unity if v falls inside the bin b and zero 
elsewhere. The characteristic error resulting from this 
simplification should be of the order of a few times 
Alogi^,, the log-spacing between bins, since the recom- 
bination timescale is a few times shorter than the Hub- 
ble tim^] In order to reduce the error induced by this 
simplification, we enforce that the bin centers coincide 
with the lowest order transition they may contain (and 
are logarithmically spaced otherwise). As can be seen 



4 Other truncation schemes could be imagined, such as assuming 
that the excited states above some threshold are in Saha equi- 
librium with the continuum; in practice, as long as a clear con- 
vergence is exhibited with increasing ri max , we need not worry 
about the detailed truncation prescription. 



5 This can be understood from Eq. jl3| : a fractional error 
A log v in the rest-frame frequency translates into the same frac- 
tional error in the redshift of emission, hence a fractional error 
A log j// (Ht) on the emissivity, where r ~ H~ l /few is the char- 
acteristic time of evolution of the populations. 
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FIG. 2. Upper panel: Weighted average of effective conductances \Q^i n + \Q 2 ^, n at T = 3800 K, as a function of transition 
energy hv n > n , for n max = 100 (open circles) and n max = 200 (dots). The nearly vertical families of points correspond to a given 
series n! — > n with fixed n, such as the Balmer series n' — > 2, the Paschen series n — > 3, etc., with n increasing from right 
to left. The families of diagonal lines correspond to a given order (n + An) — > n with fixed An, such as the a— transitions 
(n + 1) — > n, the /3— transitions (n + 2) — > n etc., with lower conductances for higher order transitions. Lower panel: Weighted 

average of effective free-bound conductances per log-frequency interval \v-^- + f v at T = 3800 K, as a function of photon 
energy hu, for n max = 100 (dotted line) and n max = 200 (solid line). The edges correspond to photoionization thresholds from 
various shells. We only show the effective conductances for these two values of n max in order not to cluster the figure; proper 
convergence with n max is discussed in Section [rv] 
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from Fig. [2j effective conductances indeed decrease with 
increasing transition order (see also Table 1 of Ref. 18). 

With this discretization, the emissivity for bound- 
bound and free-bound transitions between excited states 
becomes 



1 . I high- n 
^ v I used 



"H 



hv 

47T 



b n HX e Xp 



AT 



Qt A Xi js(v - v b ).(43) 



i=2s,2p 



the effective recombination coefficients Ai(T m ,T T ), pho- 
toionization rates Z?i(T r ) and transition rates lZi^j(T r ) 
for i,j = 2s, 2p. All the complication of the recombi- 
nation computation then resides in the slow transitions 
to the ground state, where a proper time-dependent ra- 
diative transfer calculation is required for percent-level 
precision. Here we are not worried about such subtleties, 
and use the following simple prescriptions for decay rates 
to the ground state. For the two-photon transitions from 
2s, we neglect stimulated decays and absorptions of non- 
thermal photons, so the net decay rate is 



C. Lyman-Q and 2s — Is emission 



I2 7 



X2s 



x ls e 



-E 21 /kT r 



(48) 



The net emissivity in the Lyman-a line is given by 



"H 



7- (%2p - 3&i s e 
x P CS cA 2pi i s 6{iy - v Lya ) 



(44) 



where P csc is the escape probability for the optically thick 
Lyman-a line. In the Sobolev approximation (see for 
example Ref. [41] for a detailed derivation) , it is given by 



8ttHv 3 



3c 3 nH_x ls A 2 ps s ' 
The 2s — Is two-photon emissivity is given by 



(45) 



3v 



2-; 



hv dA 2s ,_ 
47r dv 



[X2s(l + f»>)(l + U)-X ls frfr, 

(46) 

where 1/ = z/L ya — v and /„, fr are the values of the pho- 
ton occupation number at v, v' . This expression properly 
accounts for stimulated decays and absorption of non- 
thermal photons (emitted in the Lya line for example). 
In order to be consistent with our simple "standard" 
treatment of two-photon decays (see next section), we 
shall neglect these two effects here and use the approxi- 
mate expression 



■7" 1 27 hv dA 



"H 



■in 



2s, Is 

dp 



%2s 



x ls e 



(47) 



We approximated the differential two-photon decay rate 
dA.2 S ,\s/dv with the fitting formula of Ref. [12] . 

Note that neither Eq. (44 1 nor Eq. ( |47| are accurate 
at the percent level. If needed, it would be relatively 
straightforward to include the appropriate corrections. 



D. Fast part of the computation 

In this section we briefly recall how the recombination 
history can be very efficiently computed with the EMLA 
method [29] and give explicit equations for the "voltages" 
that source the emissivities. 

The EMLA formulation of the problem collapses all 
the fast and thermally-mediated interior transitions into 



where A 2s .i s ~ 8.22 s _1 is the spontaneous two-photon 
decay rate. For the net decay rate in the Lyman-a line, 
we use the Sobolev approximation, 



Xls 



Lya 



~X 2 p 



Lya 



Rhya 



x 2p - 3xi s e 



-E 21 /kT r 



where Rh ya = ^bp.is^csc is the net decay rate in the 
Ly-a line, accounting for the small escape probability of 
resonant photon, given in Eq. (45 1. 



The steady-state rate equations for 2s and 2p then read 
~ x 2s = nnx e XpA 2 s + xi s e~ E2l/kTl A 2S: is 

+ X2p'R-2p^2s ~ ^2sX2s, 

~ x 2s = n^x e x p A 2 p + 3x ls e~ E2l/kTl R Lya 

+ X2sT^2s-f2p ~ ^2pX2pj 



(50) 



(51) 



where the effective inverse lifetimes of the 2s and 2p 
states are given by 



^2s — &2s 
T2p = B2p 



7Z-2s->2p 
7^2p->2s 



A2s,ls, 
Rhya- 



(52) 
(53) 



This simple 2 by 2 system (which is exact in the limit 
that 2s and 2p are the only interface states) can be solved 
analytically. Using detailed balance relations satisfied by 
the effective rates, we find that the departures from Saha 
equilibrium are given by 



Ax 2s 



Ax 2p 



S2s 



K 2p - 



r 2s - R 

S2p + S 2 



2s-i-2p~ 

n 2 



TZ 2p 



■+2p 



-p rj~f *R- 2 s^ 2 p 

1 2p — l<-2p^,2s — fr 



where we have defined 

s 2s = n H x e x p AA 2s + Ax ls e~ E2l/kT 'A 2sAs , 

S2p = n-HXeXp AA 2 p + 3 Axis e" £2l/feTr i?Lya, 



(54) 
(55) 



(56) 
(57) 



where AAi = Ai(T m , T r )-Ai(T r , T r ) and Aa:i s is the de- 
parture of the ground state population from Saha equilib- 
rium with the plasma, defined as in Eq. Q. Similar ex- 
pressions can easily be obtained for the departures from 
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Boltzmann equilibrium with the ground state, needed for 
the Lyman-a and 2s — Is emissivities. 

The rate of change of the free-electron fraction can be 
obtained, for example, from the departures from Saha 
equilibrium computed above: 



[nnXeXpAAi 

=2s,2p 



(58) 



Finally, the matter temperature is obtained from the 
Compton-heating equation, 



-2HT„ 



&x e (iTa v Tf 



3(1 + x e + fn e )m e c 



(T z -T m ), (59) 



where <tt is the Thomson cross-section, a T is the radiation 
constant, m e is the electron mass and fn e is the He:H 
abundance ratio. At high redshifts where T m is locked to 
T r , one can use the quasi-steady-state solution [37j 



AT 



3if(l 



/ He )m e c 



8x e (TTa r T r 4 



(60) 



As an aside, we point out that it is straightforward 
to include the effect of dark matter annihilations in this 
system of equations (see for example Refs. [2HH3]). O nc 
simply has to add an additional photoionization term to 
the free-electron fraction evolution equation, a heating 
term to the matter-temperature evolution equation, and 
properly account for the additional excitations when solv- 
ing for the populations of the excited states. 



E. Numerical solution of the radiative transfer 
equation 

We first solve for the ionization history and matter 



temperature as described in Section IV D using the code 
HyRec in EMLA mode, which we have adapted to also 
extract the populations of the excited states (more pre- 
cisely, their departures form equilibrium) . We store these 
quantities on a fine redshift grid ranging from z = 2500 
to z = 400 for future interpolation. 

We follow the spectral distortion from z = 2500 to 
z = 400 on a constant energy range 10 _5 eV < E < 10.2 
eV, where the upper limit corresponds to the Lyman-a 
frequency. We assume a purely thermal spectrum blue- 
ward of Ly-a. Below z — 400, we set the emissivity to 
zero and freely redshift the distortion down to z = 0. 

We discretize the 2s — Is emissivity on the same grid 
as the high-n transitions emissivity. Our total discretized 
emissivity therefore has the form 



3v\ 



used 



J b S (» 



Vb) 



(61) 



At each time step, we first redshift the pre-existing 
spectral distortion from one bin to the next lower one, 
using Eq. ( 12 ). We then update it by adding the emission 
from the "lines" at frequencies i>b, as in Eq. (14). This 



procedure was used (with additional complications due to 
large optical depths and frequency diffusion) in Refs. [35l 
I57| . This method forces the timestep A log a to be no 
greater than the smallest bin separation A log v. 

In our fiducial computation, we have grouped effective 
conductances in bins of width A log v = 10~ 2 (and we re- 
call that the central frequency assigned to each bin is cho- 
sen to coincide with the frequency of the lowest-order line 
it contains), and used a time-step A log a = 5 x 10~ 3 . We 
have checked that reducing the bin width and timestep 
by a factor of 10 leads to maximum changes of at most 
a percent over the whole range of frequencies considered, 
with a root-mean-square difference of the order of 0.15%. 
We also checked that the spectrum is converged with re- 
spect to the redshift range over which it is computed - 
this stems from the fact that the emissivities are rela- 
tively sharply peaked around z ~ 1300, as we show in 
the next section. 



F. Results 

All quantities shown in this section are evaluated for a 
flat universe with fiducial cosmological parameters con- 
sistent with the latest WMAP results pQ: T cmb = 2.728 

k, n b h 2 = 0.022, n m h 2 = o.i3, n A h 2 = 0.34, y Hc = 

0.24, A^cff = 3.04. 

In Fig. [3j we show the number of photons emitted 
per hydrogen atom per logarithmic redshift interval or 
per logarithmic interval of observed frequency, for the 
first few transitions of the Balmer series, and for the a- 
transitions of the first few series. We see that in general 
the number of emitted photons decreases rapidly with the 
order of the transition within a series, and decreases as 
well (but less rapidly) for higher series. Note that some 
transitions may show absorption, as the H/3 transition 
P) . We find that a total of 0.63, 0.019, 0.036, 0.32 and 
0.13 photons are emitted per hydrogen atom in the Ha, 
H/3, H7, Pa and Bra transitions, respectively. 

Figure [4] shows the convergence of the high-n bound- 
bound and free-bound spectra with n max . We find that 
the fractional difference in the total spectrum between 



250 and n n 



500 is less than a percent for 



v > 0.5 GHz. Nothing formally limits us from going 
beyond n max = 500 with our method, since the tabu- 
lation of effective conductances needs to be done only 
once. However, consistently computing the spectrum at 
low frequencies would also require accounting for free- 
free absorption and collisional transitions [TSJ [25] , which 
we do not include here. We therefore limit ourselves to 
"max = 500 for now, keeping in mind that the spectrum 
obtained is not accurate below a few tenths of GHz. 

Figure [5] shows the total recombination spectrum to- 
day, as well as its subcomponents: bound-bound, free- 
bound, two-photon emission from 2s — Is decays, and 
Lyman-a emission. Note that in the present work we 
have not accounted for any of the radiative transfer ef- 
fects recently investigated with the purpose of obtaining 
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FIG. 3. Number of photons per hydrogen atom emitted per logarithmic redshift interval (or equivalently per logarithmic 
interval of today's observation frequency), for several bound-bound transitions. In the left panel, this quantity is plotted as a 
function of emission redshift. In the right panel, it is plotted as a function of observed frequency today. The two are related 

by U obs = fem/(l + *em). 



high-accuracy recombination histories (see for example 
Refs. [37J HU |4"4TI4"5] and many more references therein) . 
It is to be expected that these corrections will lead to a 
few percent corrections to the recombination spectrum; 
nevertheless, we are far from even detecting the recombi- 
nation spectrum, and such refinements are not yet needed 
for this purpose. The effective conductance method is, 
moreover, oblivious to all the complications that may 
occur between interface states and the ground state, and 
could be very easily adapted to include these effects^] 
If needed for some reason, a high-accuracy Lyman-lines 
and 2s — Is spectrum can be extracted from the modern 
recombination codes HyRec [35] and CosmoRecJ^] [50] . 
that do account for these radiative transfer effects. 

Let us point out that even though we have performed 
all computations with the correct matter temperature, 
we found that setting T m = T r (and doing so consis- 
tently, including when computing the ionization history 
and departures from Saha equilibrium), leads to an error 
of at most 0.4% for v > 0.1 GHz. It would therefore 
be sufficient, at the percent level of accuracy, to assume 
T m = T r for the spectrum computation (this assumption 
is not valid if one wishes to compute the low-redshift tail 
of the recombination history, when the two temperatures 
may differ significantly). 

Finally, we have compared our results with those of 
Ref. PH] and found a very good agreement. More thor- 
ough comparisons will be made once we implement emis- 
sion from helium as well. 



The recurring time for computing a full recombina- 
tion spectrum, independent of n max used for the effective 
conductances, is approximately 0.1 second on a standard 
workstation. This is four to six orders of magnitude faster 
than the standard multilevel atom methocr] depending 
on the value used for n max [25] . 

V. CONCLUSIONS AND FUTURE 
DIRECTIONS 

In this paper we have introduced a new and highly ef- 
ficient method to compute the primordial recombination 
spectrum. Our method relies on the factorization of the 
problem in a computationally expensive but cosmology- 
independent part and a very fast cosmology-dependent 
part - the computational efficiency of the latter part is 
itself due to a similar method of factorization, the ef- 
fective multilevel atom method, introduced in an earlier 
work [22] • The recurring cost of a spectrum computation 
with this method is a fraction of a second, at least four 
orders of magnitude faster than computations carried out 
with the standard multilevel atom method. 

The computations carried out here relied on several 
simplifying assumptions. First, we have restricted our- 
selves to the "standard" (i.e. simplified) transitions from 
2s and 2p to the ground state, and have neglected a whole 
suite of radiative transfer effects that were shown to be 



A subtlety might arise when dealing with two-photon decays 
from higher levels, as one must avoid double-counting of the low- 
frequency photons. 

http: / / www.cita.utoronto.ca/ ~jchluba 



One can speed up the computation of the spectrum by a factor of 
~ 10 if precomputing the recombination history with the EMLA 
method 29 and then computing the spectrum with the standard 
MLA method but with a larger timestep (J. Chluba, private 
communication). This would still be a few orders of magnitude 
slower than the method we present here. 
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FIG. 4. Bound-bound (not including 2s — > Is and Ly-a photons) and free-bound spectra for various values of the cutoff 
principal quantum number n max . 
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FIG. 5. Total spectral distortion created by hydrogen emission (neglecting the influence of helium), as well as individual 
contributions from various processes, using n max = 500. The "bound-bound" curve only accounts for n' — > n > 2 transitions. 
We display A/„ in the same units as Refs. [12] [H] for an easier comparison by eye. 
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important for a highly accurate recombination history 
and are efficiently computed by modern recombination 
codes [35J ED] ■ To our knowledge, no computation of the 
primordial recombination spectrum thus far has also in- 
cluded all these effects consistently. Including them with 
our formalism should not raise any major issues, but they 
are expected to induce corrections to the spectrum at the 
level of a few percent at most (though in principle one 
should check this statement precisely). The second as- 
pect we neglected is collisional transitions. They were 
shown to modify the spectrum mostly in the low fre- 
quency regime v < 0.1 GHz, where the recombination 
distortions are very smooth and therefore probably dif- 
ficult to observe [28]. Finally, we neglected broadening 
of the lines by electron scattering, as well as free-free ab- 
sorption. The former effect should increase the width of 
the broad spectral features by less than a percent of their 
central frequency, and the latter strongly affects the spec- 
trum at frequencies below ~ 0.1 GHz |19) . Since all these 
effects are small in the frequency regime > 0.1 GHz, we 
do not consider them a priority before the observability 
of the recombination spectrum is solidly established. 

In this work we have only dealt with the line and con- 
tinuum emission from recombining hydrogen atoms, and 
have neglected the effect of helium on the recombina- 
tion spectrum (but we did account for helium when com- 
puting the ionization history of the plasma). The pres- 
ence of about 0.08 helium nuclei per hydrogen nucleus 
will modify our results in the two following ways. First, 
photons emitted in the Hel 2 X P— 1 X S line (the equiva- 
lent of the hydrogen Ly-a line) with an energy of 21.2 
eV can photoionize some of the few neutral hydrogen 
atoms already present at the epoch of HcII^ Hel recom- 
bination [5lTf54] . This effect is already included in our 
code when computing the ionization fraction of helium 
[35j , but we have not accounted here for the small depar- 
ture from Saha equilibrium that it induces for hydrogen, 
and the resulting line emission. This process was studied 
in detail in Ref. [53], who showed that it leads to pre- 
recombination emission features in the hydrogen spec- 
trum. Inclusion of this effect with our method presents 
no additional complication: one only needs to properly 
compute the small departures from equilibrium in hy- 
drogen, given the helium ionizing photons. The second 
and more direct consequence of having helium nuclei is 
that they too recombined, and in the process, emitted a 
few nonthermal photons per helium nucleus. Because he- 
lium recombined in two stages (at redshifts z ~ 6000 for 
HellF- > Hell recombination and z ~ 2000 for Hell— s- Hel 
recombination), this results in ~ 0.16 helium recombina- 
tion photons per hydrogen recombination photon, which 



is a non-negligible contribution to the spectrum, and can 
be used to probe the primordial helium abundance be- 
fore the first stars formed [50]. Computing the emission 
from helium recombination can be achieved with the ex- 
act same method, the only additional work would be to 
compute effective conductances for Hel (those for Hell 
can be easily obtained from the ones computed for hy- 
drogen by simple rescalings). Such additions are essential 
and we shall include them in the near future. 

An additional aspect that we have not treated here is 
the processing of pre-existing spectral distortions by the 
recombining atoms. These can lead to large enhance- 
ments of the line emission, even if the initial smooth 
spectral distortion is tiny [23]. Formally, one can gener- 
alize our method to account for non-thermally mediated 
transitions, as long as the underlying spectrum can be 
described by a small set of parameters. We defer such 
a study to future work. Once the latter two aspects are 
implemented, we shall release a public code to efficiently 
compute the recombination spectrum. 

Besides having a certain aesthetic appeal, we believe 
that the factorization method presented here will prove 
very useful in order to quantitatively study the infor- 
mation content of the recombination spectrum. Most of 
our knowledge about the early universe currently comes 
from the study of spatial variations of the thermal photon 
background, through CMB anisotropy observations. An 
exciting avenue is to probe the spectral variations of the 
CMB, which might inform us about the universe before 
the last-scattering redshift. A fast and accurate code to 
compute the primordial spectrum will be of great help to 
study precisely if, what and how new information can be 
gained from the primordial spectrum. 
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